/**
 * @file MPI/cpi.cpp
 *
 * @license
 * Copyright (c) 2009-2018
 * Fraunhofer Institute for Algorithms and Scientific Computing SCAI
 * for Fraunhofer-Gesellschaft
 *
 * This file is part of the SCAI framework LAMA.
 *
 * LAMA is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.
 *
 * LAMA is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
 * more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with LAMA. If not, see <http://www.gnu.org/licenses/>.
 * @endlicense
 *
 * @brief Simple mpi example calculating pi.
 * @author Jan Ecker
 * @date 20.03.2013
 */

#include "mpi.h"
#include <cstdio>
#include <cmath>

static double f(double);

double f( double a )
{
    return (4.0 / (1.0 + a*a));
}

int main( int argc,char *argv[] )
{
    int done = 0, n, myid, numprocs, i;
    double PI25DT = 3.141592653589793238462643;
    double mypi, pi, h, sum, x;
    double startwtime = 0.0, endwtime;
    int  namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);
    MPI_Get_processor_name(processor_name,&namelen);

    fprintf(stdout,"Process %d of %d on %s\n", myid, numprocs, processor_name);

    done = 10;
    n = 0;
    while ( !done )
    {
        if (myid == 0)
        {
            /*
            printf("Enter the number of intervals: (0 quits) ");
            scanf("%d",&n);
            */

	        if (n==0) n=10000; else n=0;

	        startwtime = MPI_Wtime();
        }

        MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
        
        if (n == 0)
        {
            done -= 1;
        }
        else
        {
            h   = 1.0 / (double) n;
            sum = 0.0;

	        /* A slightly better approach starts from large i and works back */
            for (i = myid + 1; i <= n; i += numprocs)
            {
                x = h * ((double)i - 0.5);
                sum += f(x);
            }
            mypi = h * sum;

            MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

            if (myid == 0)
	        {
                printf("pi is approximately %.16f, Error is %.16f\n", pi, fabs(pi - PI25DT));
		        endwtime = MPI_Wtime();
		        printf("wall clock time = %f\n", endwtime-startwtime);	       
		        fflush( stdout );
	        }
        }
    }
    MPI_Finalize();
    return 0;
}
